B10-L10-int and B10-P10-int DE (neurons)

library(modplots)

library(tidyverse)

── Attaching core tidyverse packages ───────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(Seurat)

Loading required package: SeuratObject
Loading required package: sp

Attaching package: ‘SeuratObject’

The following objects are masked from ‘package:base’:

    intersect, t

library(plotly)

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

# Mitochondrial genes
mt <- c("ENSGALG00000035334", #COX3
        "ENSGALG00000032456", #COII
        "ENSGALG00000032142", #MT-CO1
        "ENSGALG00000032079", #MT-CYB
        "ENSGALG00000037838", #ND6
        "ENSGALG00000029500", #ND5
        "ENSGALG00000036229", #MT-ND4
        "ENSGALG00000042478", #ND4L
        "ENSGALG00000030436", #ND3
        "ENSGALG00000041091", #MT-ATP6
        "ENSGALG00000032465", #MT-ATP8
        "ENSGALG00000043768", #MT-ND2
        "ENSGALG00000042750") #MT-ND1

# volcanoplot thresholds
p.adj <- 0.05
l2fc <- 0.5

B10-L10-int

# seurat object
my.se <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_int_seurat_250723.rds")
# cluster labels from B10int and L10int
ctrl_lumb_int_combined_labels <- readRDS("~/spinal_cord_paper/annotations/ctrl_lumb_int_combined_labels.rds")

identical(colnames(my.se), rownames(ctrl_lumb_int_combined_labels))

[1] TRUE

my.se$annot_sample  <- ctrl_lumb_int_combined_labels$annot_sample

my.se@active.assay <- "RNA"

intra cluster DE analysis

We do DE analysis per cluster, contrasting B10 and L10 samples:

markers <- list()
numbers <- list()
composition <-list()

for (i in seq(levels(Idents(my.se)))) {
  # subset for individual clusters
  mn.se <- subset(x = my.se, idents = levels(Idents(my.se))[i])
  mn.se$sample <- str_extract(mn.se$orig.ident, "ctrl|lumb")
  
  composition[[i]] <- mn.se[[]] %>% 
    select(sample, annot_sample)
  
  Idents(mn.se) <- "sample"
  
  tmp_markers <- FindMarkers(mn.se,
                             ident.1 = "ctrl",
                             only.pos = FALSE, 
                             min.pct = 0.25, 
                             logfc.threshold =  0.2,
                             latent.vars = c("CC.Difference.seurat"), 
                             test.use = "MAST", 
                             assay = "RNA")
  # cell numbers per sample
  numbers[[i]] <- data.frame(table(mn.se$sample))
  
  tmp_markers <- tmp_markers %>%
    rownames_to_column("Gene.stable.ID") %>% 
    left_join(gnames)
  
  markers[[i]] <- tmp_markers
}

names(markers) <- paste0("cl-", levels(Idents(my.se)))
names(numbers) <- paste0("cl-", levels(Idents(my.se)))
names(composition) <- paste0("cl-", levels(Idents(my.se)))

# bind lists into data frames
lumb_markers <- bind_rows(markers, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))
lumb_numbers <- bind_rows(numbers, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))
lumb_composition <- bind_rows(composition, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))

saveRDS(lumb_markers, "~/spinal_cord_paper/data/Gg_ctrl_lumb_int_markers.rds")
saveRDS(lumb_numbers, "~/spinal_cord_paper/data/Gg_ctrl_lumb_int_numbers.rds")
saveRDS(lumb_composition, "~/spinal_cord_paper/data/Gg_ctrl_lumb_int_composition.rds")

lumb_markers <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_int_markers.rds") %>% 
  mutate(clust_id = str_remove(cluster, "^cl-"))

Error in readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_int_markers.rds") %>%  : 
  could not find function "%>%"

Plot the cluster compositions and the number of marker genes.

DimPlot(my.se, label = TRUE, reduction = "tsne")

ggplot(data = lumb_numbers, aes(x = Var1, y = Freq, group = cluster, label = Freq)) +
  geom_col() +
  facet_wrap("cluster", nrow = 3) +
  geom_text(nudge_y = 50, size = 3) +
  ggtitle("Cluster composition by sample (nCells)")

ggplot(data = lumb_markers %>% filter(p_val_adj < 0.05 & avg_log2FC > 0.5), aes(x = cluster, group = cluster)) +
  geom_bar() +
  facet_wrap("cluster", nrow = 3,scales = "free_x") +
  ggtitle("Number of sig. DE genes")

NA
NA
NA

DA clusters from miloR

We select the clusters with the highest amount of neighbourhoods with DA of cells. Those are clusters 3 (neurons), 21, 22 (Glycogen body), and 30 motor neurons.

# vector of clusters
select_clusters <- c(3,21,22,30)

select_lumb_markers <- filter(lumb_markers, clust_id %in% select_clusters)
select_lumb_numbers <- filter(lumb_numbers, clust_id %in% select_clusters)
select_lumb_composition <- filter(lumb_composition, clust_id %in% select_clusters)

Barplots

Barplots showing number of cells from B10 or L10, and the contributions from the individual B10int and L10int clusters:

ggplot(data = select_lumb_numbers, aes(x = Var1, y = Freq, group = cluster, label = Freq)) +
  geom_col() +
  facet_wrap("cluster", nrow = 1) +
  geom_text(nudge_y = 10) +
  ggtitle("Cluster composition by sample (nCells)")

toplot <- select_lumb_composition %>% 
  mutate(annot_sample = str_extract(annot_sample, "\\d{1,2}_.{1}(?=trl|umb)")) %>% 
  mutate(annot_sample = factor(annot_sample)) %>% 
  group_by(sample, cluster) %>%
  count(annot_sample) %>% 
  mutate(label_ypos = cumsum(n)- 0.5*n)

BL_select_barplot <- ggplot(data=toplot, aes(x=sample, y=n, fill=fct_rev(annot_sample))) +
  geom_bar(stat="identity", color = "black") +
  geom_text(aes(y=label_ypos, label=annot_sample), 
            color="black", size=3.5) +
  facet_wrap("cluster", nrow = 1) +
    theme(legend.position="none") +
  ggtitle("BL10int cluster contributions")

BL_select_barplot

pdf("~/spinal_cord_paper/figures/Fig_4_BL10int_high_DA_clust_contribution.pdf")

BL_select_barplot
dev.off()

png 
  2 

Volcanoplots

toplot <- select_lumb_markers %>% 
  mutate(gene_type = case_when(
    avg_log2FC >= 0.5 & p_val_adj <= 0.05 ~ "ctrl",
    avg_log2FC <= -0.5 & p_val_adj <= 0.05 ~ "lumb",
    TRUE ~ "ns")
  )

cols <- c(ctrl = "black",
          lumb = "#419c73",
          ns = "grey")

shapes <- c(ctrl = 21,
            lumb = 21,
            ns = 20)

volplot <- ggplot(data = toplot,
                  aes(x = avg_log2FC,
                      y = -log10(p_val_adj),
                      label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                  )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  facet_wrap("cluster", nrow = 1) +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes")

volplot

# ggplotly(volplot, width = 1000, height = 600)

Plots without HOX and mitochondrial genes:

# filter out HOX and MT genes
toplot_nomt <- toplot %>% 
  filter(!grepl("^HOX", Gene.name)) %>% 
  filter(!Gene.stable.ID %in% mt)

volplot_nomt <- ggplot(data = toplot_nomt,
                  aes(x = avg_log2FC,
                      y = -log10(p_val_adj),
                      label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                  )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  facet_wrap("cluster", nrow = 1) +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes (without MT and HOX genes")

volplot_nomt

# ggplotly(volplot_nomt, width = 1000, height = 600)

volplot_nomt +
    ggrepel::geom_text_repel(aes(label = ifelse(gene_type == 'ns',
                                                NA,
                                                Gene.name)),
                             size = 3,
                             color = "black")

Warning: Removed 8233 rows containing missing values or values outside the scale range (`geom_text_repel()`).
Warning: ggrepel: 6 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 208 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 534 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 52 unlabeled data points (too many overlaps). Consider increasing max.overlaps

No MN DE genes

Cluster 30 (Consisting of 161 B10 vs 13 L10 cells) shows no MN related markers. Seemingly, those 13 cells are indeed motor neurons. This is supported by the fact of their expression of TUBB3, FOXP1, and SLC18A3.

# Motor neuron cluster
mn.se <- subset(x = my.se, idents = 30)
mn.se$sample <- str_extract(mn.se$orig.ident, "ctrl|lumb")

mn.se@active.assay <- "integrated"

mn_markers <- gnames[grepl("TUBB3|FOXP1|SLC18A3", gnames$Gene.name),]
mn_markers

VlnPlot(mn.se, group.by = "sample", mn_markers$Gene.stable.ID, cols = c("darkgrey", "#419c73"))

Low DA clusters from miloR

In contrast to the clusters above, we now select clusters 17, 18, and 19. The show a very low amount of DA neighbourhoods, are of similar size and represent neuronal, MFOL, and NPC populations.

# vector of MFOL clusters
low_clusters <- c(17, 18, 19)

low_lumb_markers <- filter(lumb_markers, clust_id %in% low_clusters)
low_lumb_numbers <- filter(lumb_numbers, clust_id %in% low_clusters)
low_lumb_composition <- filter(lumb_composition, clust_id %in% low_clusters)

Barplots

Barplots showing number of cells from B10 or L10, and the contributions from the individual B10int and L10int clusters:

ggplot(data = low_lumb_numbers, aes(x = Var1, y = Freq, group = cluster, label = Freq)) +
  geom_col() +
  facet_wrap("cluster", nrow = 1) +
  geom_text(nudge_y = 10) +
  ggtitle("Cluster composition by sample (nCells)")

toplot <- low_lumb_composition %>% 
  mutate(annot_sample = str_extract(annot_sample, "\\d{1,2}_.{1}(?=trl|umb)")) %>% 
  mutate(annot_sample = factor(annot_sample)) %>% 
  group_by(sample, cluster) %>%
  count(annot_sample) %>% 
  mutate(label_ypos = cumsum(n)- 0.5*n)

BL_low_barplot <- ggplot(data=toplot, aes(x=sample, y=n, fill=fct_rev(annot_sample))) +
  geom_bar(stat="identity", color = "black") +
  geom_text(aes(y=label_ypos, label=annot_sample), 
            color="black", size=3.5) +
  facet_wrap("cluster", nrow = 1) +
    theme(legend.position="none") +
  ggtitle("From which clusters do cells come?")

BL_low_barplot

pdf("~/spinal_cord_paper/figures/Fig_4_BL10int_low_DA_clust_contribution.pdf")

BL_low_barplot
dev.off()

png 
  2 

Volcanoplots

toplot <- low_lumb_markers %>% 
  mutate(gene_type = case_when(
    avg_log2FC >= 0.5 & p_val_adj <= 0.05 ~ "ctrl",
    avg_log2FC <= -0.5 & p_val_adj <= 0.05 ~ "lumb",
    TRUE ~ "ns")
  )

cols <- c(ctrl = "black",
          lumb = "#419c73",
          ns = "grey")

shapes <- c(ctrl = 21,
            lumb = 21,
            ns = 20)

volplot <- ggplot(data = toplot,
                  aes(x = avg_log2FC,
                      y = -log10(p_val_adj),
                      label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                  )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  xlim(c(-3.5, 3.5)) +
  facet_wrap("cluster", nrow = 1) +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes")

volplot

# ggplotly(volplot, width = 800, height = 500)

Plots without HOX and mitochondrial genes:

# filter out HOX and MT genes
toplot_nomt <- toplot %>% 
  filter(!grepl("^HOX", Gene.name)) %>% 
  filter(!Gene.stable.ID %in% mt)

volplot_nomt <- ggplot(data = toplot_nomt,
                  aes(x = avg_log2FC,
                      y = -log10(p_val_adj),
                      label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                  )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  xlim(c(-1.5, 1.5)) +
  facet_wrap("cluster", nrow = 1, scales = "free_y") +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes (without MT and HOX genes)")

volplot_nomt

# ggplotly(volplot_nomt, width = 800, height = 500)

p1 <- volplot +
    ggrepel::geom_text_repel(aes(label = ifelse(gene_type == 'ns',
                                                NA,
                                                Gene.name)),
                             size = 3,
                             color = "black")

p2 <- volplot_nomt +
    ggrepel::geom_text_repel(aes(label = ifelse(gene_type == 'ns',
                                                NA,
                                                Gene.name)),
                             size = 3,
                             color = "black")

p1

Warning: Removed 1289 rows containing missing values or values outside the scale range (`geom_text_repel()`).
Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 15 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 34 unlabeled data points (too many overlaps). Consider increasing max.overlaps

p2

Warning: Removed 1264 rows containing missing values or values outside the scale range (`geom_text_repel()`).

pdf("~/spinal_cord_paper/figures/Fig_4_BL10int_low_DA_volplots.pdf", width = 8, height = 4)

p1

Warning: Removed 1289 rows containing missing values or values outside the scale range (`geom_text_repel()`).
Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 16 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 35 unlabeled data points (too many overlaps). Consider increasing max.overlaps

p2

Warning: Removed 1264 rows containing missing values or values outside the scale range (`geom_text_repel()`).
Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider increasing max.overlaps

dev.off()

png 
  2 

B10-P10-int

# seurat object
my.se <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_poly_int_seurat_250723.rds")
# cluster labels from B10int and L10int
ctrl_poly_int_combined_labels <- readRDS("~/spinal_cord_paper/annotations/ctrl_poly_int_combined_labels.rds")

identical(colnames(my.se), rownames(ctrl_poly_int_combined_labels))

[1] TRUE

my.se$annot_sample  <- ctrl_poly_int_combined_labels$annot_sample

my.se@active.assay <- "RNA"

intra cluster DE analysis

We do DE analysis per cluster, contrasting B10 and P10 samples:

markers <- list()
numbers <- list()
composition <- list()

for (i in seq(levels(Idents(my.se)))) {
  # subset for individual clusters
  mn.se <- subset(x = my.se, idents = levels(Idents(my.se))[i])
  mn.se$sample <- str_extract(mn.se$orig.ident, "ctrl|poly")
  
  composition[[i]] <- mn.se[[]] %>% 
    select(sample, annot_sample)
  
  Idents(mn.se) <- "sample"
  
  tmp_markers <- FindMarkers(mn.se,
                             ident.1 = "ctrl",
                             only.pos = FALSE, 
                             min.pct = 0.25, 
                             logfc.threshold =  0.2,
                             latent.vars = c("CC.Difference.seurat"), 
                             test.use = "MAST", 
                             assay = "RNA")
  # cell numbers per sample
  numbers[[i]] <- data.frame(table(mn.se$sample))
  
  tmp_markers <- tmp_markers %>%
    rownames_to_column("Gene.stable.ID") %>% 
    left_join(gnames)
  
  
  markers[[i]] <- tmp_markers
}

names(markers) <- paste0("cl-", levels(Idents(my.se)))
names(numbers) <- paste0("cl-", levels(Idents(my.se)))
names(composition) <- paste0("cl-", levels(Idents(my.se)))

# bind lists into data frames
poly_markers <- bind_rows(markers, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))
poly_numbers <- bind_rows(numbers, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))
poly_composition <- bind_rows(composition, .id = "cluster") %>% 
  mutate(cluster = factor(cluster, levels = paste0("cl-", levels(Idents(my.se)))))

saveRDS(poly_markers, "~/spinal_cord_paper/data/Gg_ctrl_poly_int_markers.rds")
saveRDS(poly_numbers, "~/spinal_cord_paper/data/Gg_ctrl_poly_int_numbers.rds")
saveRDS(poly_composition, "~/spinal_cord_paper/data/Gg_ctrl_poly_int_composition.rds")

# load the DE data
poly_markers <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_poly_int_markers.rds") %>% 
  mutate(clust_id = str_remove(cluster, "^cl-"))
poly_numbers <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_poly_int_numbers.rds") %>% 
  mutate(clust_id = str_remove(cluster, "^cl-"))
poly_composition <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_poly_int_composition.rds") %>% 
  mutate(clust_id = str_remove(cluster, "^cl-"))

Plot the cluster compositions and the number of marker genes.

DimPlot(my.se, label = TRUE, reduction = "tsne")

ggplot(data = poly_numbers, aes(x = Var1, y = Freq, group = cluster, label = Freq)) +
  geom_col() +
  facet_wrap("cluster", nrow = 3) +
  geom_text(nudge_y = 50, size = 3) +
  ggtitle("Cluster composition by sample (nCells)")

ggplot(data = poly_markers %>% filter(p_val_adj < 0.05 & avg_log2FC > 0.5), aes(x = cluster, group = cluster)) +
  geom_bar() +
  facet_wrap("cluster", nrow = 3,scales = "free_x") +
  ggtitle("Number of sig. DE genes")

Fig5 candidates

# vector of neuronal clusters
fig5_clusters <- c(5,11,18,23,25)

fig5_poly_markers <- filter(poly_markers, clust_id %in% fig5_clusters)
fig5_poly_numbers <- filter(poly_numbers, clust_id %in% fig5_clusters)
fig5_poly_composition <- filter(poly_composition, clust_id %in% fig5_clusters)

Barplots

Barplots showing number of cells from B10 or P10, and the contributions from the individual B10int and P10int clusters:

ggplot(data = fig5_poly_numbers, aes(x = Var1, y = Freq, group = cluster, label = Freq)) +
  geom_col() +
  facet_wrap("cluster", nrow = 1) +
  geom_text(nudge_y = 10) +
  ggtitle("Cluster composition by sample (nCells)")

toplot <- fig5_poly_composition %>% 
  mutate(annot_sample = str_extract(annot_sample, "\\d{1,2}_.{1}(?=trl|oly)")) %>% 
  mutate(annot_sample = factor(annot_sample)) %>% 
  group_by(sample, cluster) %>%
  count(annot_sample) %>% 
  mutate(label_ypos = cumsum(n)- 0.5*n)

BP_neuron_barplot <- ggplot(data=toplot, aes(x=sample, y=n, fill=fct_rev(annot_sample))) +
  geom_bar(stat="identity", color = "black") +
  geom_text(aes(y=label_ypos, label=annot_sample), 
            color="black", size=3.5) +
  facet_wrap("cluster", nrow = 1) + 
    theme(legend.position="none") +
  ggtitle("From which clusters do cells come?")

BP_neuron_barplot 

Volcanoplots

toplot <- fig5_poly_markers %>% 
  mutate(gene_type = case_when(
    avg_log2FC >= 0.5 & p_val_adj <= 0.05 ~ "ctrl",
    avg_log2FC <= -0.5 & p_val_adj <= 0.05 ~ "poly",
    TRUE ~ "ns")
  )

cols <- c(ctrl = "black",
          poly = "goldenrod3",
          ns = "grey")

shapes <- c(ctrl = 21,
            poly = 21,
            ns = 20)


volplot <- ggplot(data = toplot,
                       aes(x = avg_log2FC,
                           y = -log10(p_val_adj),
                           label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                       )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  facet_wrap("cluster", nrow = 2, scales = "free") +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes")

ggplotly(volplot, width = 1000, height = 600)

Warning in dev_fun(filename = tempfile(), width = width %||% 640, height = height %||%  :
  unable to open connection to X11 display ''
Error in .External2(C_X11, paste0("png::", filename), g$width, g$height,  : 
  unable to start device PNG

Plots without mitochondrial genes:

toplot_nomt <- toplot %>% 
  filter(!Gene.stable.ID %in% mt)

volplot_nomt <- ggplot(data = toplot_nomt,
                       aes(x = avg_log2FC,
                           y = -log10(p_val_adj),
                           label = Gene.name,
                      color = gene_type,
                      shape = gene_type
                       )) +
  geom_point() +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = cols) +
  scale_shape_manual(values = shapes) +
  facet_wrap("cluster", nrow = 1) +
  ylab("-log10(padj)") +
  theme_bw()+
  ggtitle("Sig. marker genes (without MT and HOX genes")

ggplotly(volplot_nomt)

Warning in dev_fun(filename = tempfile(), width = width %||% 640, height = height %||%  :
  unable to open connection to X11 display ''
Error in .External2(C_X11, paste0("png::", filename), g$width, g$height,  : 
  unable to start device PNG

pdf("~/spinal_cord_paper/figures/Fig_5_BPoly10int_volplots.pdf", width = 12, height = 4)
pfig5

Warning: Removed 2860 rows containing missing values or values outside the scale range
(`geom_text_repel()`).
Warning: ggrepel: 151 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 6 unlabeled data points (too many overlaps). Consider increasing max.overlaps

dev.off()

null device 
          1 

# Date and time of Rendering
Sys.time()

sessionInfo()